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ABSTRACT 

The  linear  stability  of  the  Batchelor  (1964)  vortex  is  investigated.  Particular  emphasis 
is  placed  on  modes  found  recently  in  a  numerical  study  by  Khorrami  (1991).  These  modes 
have  a  number  of  features  very  distinct  from  those  found  previously  for  this  vortex,  in¬ 
cluding  (i)  exhibiting  small  growth  rates  at  large  Reynolds  numbers  and  (ii)  susceptibility 
to  destabilisation  by  viscosity.  In  this  paper  these  modes  are  described  using  asymptotic 
techniques,  producing  results  which  compare  very  favourably  with  fully  numerical  results  at 
large  Reynolds  numbers. 
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1.  Introduction 


Stability  analysis  of  streamwise  vortices  plays  an  important  role  in 
such  diverse  areas  as  wake-hazard  reduction,  combustor  optimization,  and 
turbulent  boundary- layer  structure.  Employing  the  Batchelor  vortex 
(Batchelor  1964)  for  the  mean  velocity  profile,  a  great  deal  of  effort 
has  been  directed  towards  understanding  the  stability  characteristics 
of  a  trailing-line  vortex;  the  numerical  works  of  Lessen,  Singh  & 

Paillet  (1974),  Lessen  &  Paillet  (1974)  and  Duck  &  Foster  (1980)  should  be 
mentioned.  Using  asymptotic  analysis,  the  findings  of  the  above  authors 
were  confirmed  by  many  investigators,  including  Stewartson  (1982), 

Leibovich  and  Stewartson  (1983),  Stewartson  &  Capell  (1985), 

Stewartson  &  Brown  (1985),  Duck  (1986),  and  Stewartson  &  Leibovich  (1987). 
These  asymptotic  studies  reveal  the  complex  nature  and  structures  of  the 
inviscid  modes  with  negative  azimuthal  wavenumbers.  Furthermore,  they 
showed  the  intricacies  and  difficulties  associated  with  the  numerical 
computations  of  these  instabilities.  However,  most  of  the  above  studies 
treated  only  inviscid  disturbances,  and  with  the  possible  exception  of  the 
work  of  Maslowe  &  Stewartson  (1982),  viscosity  was  believed  to  have  a 
stabilizing  influence. 

Recently,  using  a  numerical  method,  Khorrami  (1991)  found  new 
viscous  modes  of  instability  for  the  Batchelor  vortex.  The  two 
reported  modes  are  for  azimuthal  wavenumbers  which  previously  were 
thought  to  be  stable.  Furthermore,  Khorrami  found  these  modes  differ 
from  the  inviscid  disturbances  studied  previously  in  two  respects.  First, 
there  are  no  higher  modes  associated  with  them,  and  second  they  have 
growth  rates  which  are  generally  orders  of  magnitude  smaller.  In 
light  of  the  above,  it  seems  quite  unlikely  for  these  new 
instabilities  to  have  structures  similar  to  the  inviscid 
perturbations  reported  by  previous  investigators.  However,  regarding  these 
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modes,  numerical  methods  are  not  the  proper  tool  for  providing 
either  scale  and  structural  information  or  a  limiting  analysis  near  the 
neutral  curves. 

This  paper  is  an  effort  to  address  these  concerns,  as  well  as 
to  provide  firmer  grounds  for  the  existence  of  the  instability  modes 
with  positive  (and  zero)  azimuthal  wavenumbers  for  the  Batchelor  vortex.  A 
combination  of  asymptotic  and  numerical  analysis  is  presented. 
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2.  Problem  Formulation 

If  (u*.v*,w*)  denote  the  dimensional  velocity  components  in  the 
radial  (r*)  azimuthal  (0)  and  axial  (x*)  directions  respectively, 
then  the  similarity  solution  of  swirling  wake  flows  at  high  Reynolds 
numbers  due  to  Batchelor  (1964)  may  be  written 


w0*  =  d-e'11) 

r 

u0*  *  u0  ■  log  ( "Di!  ]  +  0£  «rt) 

8vx*  1  v  1  8vx* 

- 

8vx* 


(2.1) 


(2.2) 


where 


_  »Qi 


4vx 


(2.3) 


v  is  the  kinematic  viscosity  of  the  (incompressible)  fluid,  L  is  a 
constant  (akin  to  a  drag  coefficient),  Cq  is  the  circulation  at 
large  radius,  and 

Q(ti)  =  e'11  (logi)  +  ei(ri)  -  0.807} 

+  2  ei(n)  -  2e  i  (2rj) , 


where 


(2.4) 


et(ri)  =  (  dC.  (2.5) 

n  s 

Batchelor  (1964)  showed  the  term  involving  Q(t|)  in  (2.2)  is 
numerically  much  smaller  in  magnitude  than  the  other  terms,  and 
consequently  will  be  neglected.  Similar  assumptions  have  been  implemented 
by  previous  studies  on  the  stability  of  this  class  of  vortical  flow,  as 
detailed  in  the  previous  section.  Following  Lessen  et  al.  (1974),  we 
scale  velocity  by 
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C02  Uqx*  LU02 


Us  =  -  log 


(2.6) 


and  length  scales  by 


r  4vx*  li 

=  l  J  • 


(2.7) 


This  leads  to  a  non-dimensional  mean-flow  profile  given  by 

u0  .  r2 
U  =  U7*  e  • 


(2.8) 


w  =  a  (l-e'r2), 


(2.9) 


where 


q  = 

r  =  r*/rs  =  >/r\. 


(2.10) 


(2.11) 


We  now  write  the  velocity  field  as  the  sum  of  the  mean  flow  together 
with  a  small  amplitude  perturbation,  viz 
u*  =  Us(U+5u) , 
v*  =  Us5v, 

w*=Us(W+5w),  (2.12) 

whilst  the  pressure  is  written  as 


where 


P*  =  pU02tn+5p] , 


Co2U0 

n  =  -T~Qi(r). 
8vx*Us2 


(2.13) 


(2.14) 


. 

Qj(r)  =  - K  2ei (r2)  -  2ei(2r2), 


(2.15) 


A  tilde  quantity  here  represents  the  perturbation  about  the  mean 
state  and  5  is  the  (small)  perturbation  amplitude.  We  now  make  the 
further  assumption  that  U  and  W  (and  indeed  also  Uq)  are  independent 
of  x*.  This  will  generally  be  an  improper  assumption,  and  is  equivalent 
to  a  parallel  flow  approximation  (which  has  been  used  as  an 
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assumption  in  numerous,  diverse,  stability  investigations  previously). 
However,  we  justify  this  step  on  the  following  grounds.  First,  one  of 
the  primary  aims  in  this  paper  is  to  develop  asymptotic  theories  to 
compare  with  previous  numerical  results,  which  were  all  based  on 
the  same  parallel  flow  approximation.  Second,  since  to  leading  order  the 
solutions  to  which  we  concern  ourselves  turn  out  to  be  inviscid  in  form,  it 
can  be  shown  that  to  first  order  the  parallel  flow  approximation  is  a 
right  and  proper  one. 

We  now  return  to  consideration  of  the  form  to  be  taken  for  the 
perturbation  quantities.  We  write 


(u.v.w.p)  =  {F(r),iG(r),H(r),P(r)}  exp{i (ax+n8-act> , 

(2.16) 

where  a  and  n  are  the  axial  and  azimuthal  wavenumbers  respectively, 

and  c  =  cr  +  i  Cj  is  the  complex  wavespeed.  It  turns  out  that  the 

problem  remains  unaltered  if  we  use 
-r2 

U  =  e  r  (2.17) 

as  the  mean  axial  velocity  distribution,  provided  we  also  replace 
”c"  by  "-c",  "W”  by  "-W"  and  "P”  by  "-P”.  The  only 
net  effect  of  this  is  on  cr,  whilst  the  important  amplification  rate 
Cj  is  totally  unaffected. 

If  we  then  substitute  (2.12),  (2.13),  and  (2.16)  into  the  equations 
of  motion,  and  consider  terms  solely  of  0(5),  we  then  obtain 


G'  +  ~  +  aF  +  n~  =  0, 


.  G"  iG*  r  i  /  n'-fi  2  1  *  1 

r  *  I  ^  t  +  * 1 ' 9  J 


(2.18) 


n2+l 


0. 


(2.19) 
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H' ' 


1 


l^T 


H*  + 


+  a2 


(2.20) 


F'  1  1 

Ke 


+ 


}) 


F 


+  i  ijy  G  +  iaP  =  0,  (2.21) 

where  a  prime  denotes  differentiation  with  respect  to  the  radial  coordinate. 
Here  the  Reynolds  number  is  defined  as 

Re  =  Hjli  ,  (2.22) 


and 

nW 

<p  =  a(U-c)  +  2*.  (2.23) 

The  boundary  conditions  that  must  be  imposed  on  this  system  are: 
at  r  =  0, 

for  n  =  0,  G(0)  =  H(0)  =  F'(0)  =  P'(0)  =  0, 

for  n  =  ±1,  G'(0)  =  G(0)±H(0)  =  F(0)  =  P(0)  =  0, 
for  n  >  1,  F(0)  =  G(0)  =  H(0)  =  P(0)  =  0,  (2.24) 

whi  1st  as  r  -♦  », 

F(r),  G(r) ,  H(r),  P(r)  ^  0.  (2.25) 

In  the  following  sections  we  study  the  above  system  in  the  limit 
as  Rc  -»  «. 


A  number  of  papers  previously  addressed  the  inviscid  limit  of  the 

system  (2.18)  -  (2.21),  (2.24),  (2.25),  in  particular  for  modes  which 
exhibit  finite  temporal  growth  rates  (acj)  in  this  limit.  Rather  than 
concern  ourselves  with  these  modes,  we  focus  on  another  family  of  modes 
found  numerically  by  Khorrami  (1991).  These  exhibit  significantly 
diminishing  growth  rates  with  a  decrease  in  viscosity.  Inspection 
of  a  number  of  these  results,  and  others  not  reported,  suggests  the 
following  two  general  characteristics  of  these  modes  as 
Re  -»  «:  (i)  cr  =  0(1)  and  (ii)  cj  =  0(Re*1).  These  trends  strongly 

suggest  we  seek  an  asymptotic  expansion  to  our  solution  of  the  form 
(F.G.H.P)  =  {Fo.Go.Hq.Po)  +  Re'^Fi.G^Hj.P!)  +  0(Re*2),  (3.1) 

c  =  c0  +  Re*1  C!  +  0(Re-2),  (3.2) 

where  we  expect  cq  to  be  real.  Substituting  (3.1),  (3.2)  into 
(2. 18)-(2.21)  and  taking  just  0(1)  terms  yields  the  following 
(inviscid)  system  of  equations 


G0’  +  +  a  F0  +  25Q  =  o. 

(3.3) 

90  G0  +  2  =  P0' , 

(3.4) 

90  F0  +  U’G0  =  -«P0. 

(3.5) 

90  »o  +  (W'4)G0  =  ’  *7^* 

(3.6) 

nW 

where  90  =  -cq)  +  — • 

(3.7) 

Indeed,  Fq  and  Hq  may  be  eliminated  between  these  equations  to 
yield  the  following  ordinary  differential  equations  as  determined 
by  Duck  and  Foster  (1980): 
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or  symbolically 

Li  { G0 , P0 }  =  0,  (3.9) 

together  wi th 


dP0  r 

IT  =  [  *0  * 


(W2r2)' 
r39  0 


]  G0 


2llW  Pn 

~2~  P°’ 
r^<P0 


(3.10) 


or  symbolically 

l2  { Gq , P0 }  =  0.  (3.11) 

The  boundary  conditions  may  be  simply  inferred  from  (2.24)  and  (2.25). 

Equations  (3.8),  (3.10)  (and  equivalent)  have  been  investigated 
by  a  nunuer  of  authors  (e.g.  Lessen  et  al  1974,  Duck  &  Foster  1980),  in 
particular  for  complex  values  of  the  wavespeed  cq.  For  this  study, 
we  carried  out  a  similar  investigation  but  sought  real  values  of  cq. 

A  Chebyshev  spectral  collocation  method  was  employed  to  perform  the 
numerical  tasks  throughout  this  study,  since  spectral  techniques  are 
well  known  for  their  accuracy  and  fast  convergence  rate.  The 
mathematical  theory  of  such  methods  is  found  in  Gottlieb  &  Orszag 
(1977)  and  Gottlieb,  Hussaini  &  Orszag  (1984)  and  is  not  presented  here. 
Its  implementation  for  the  stability  of  swirling  flows  is  given  in  detail 
by  Khorrami,  Malik  &  Ash  (1989),  and  readers  are  referred  to  that 
paper  for  further  information.  Briefly,  the  method  consists  of 
expanding  each  perturbation  eigenfunction  in  a  truncated  Chebyshev 
series,  for  example 


N 

G(5)  =  I  akTk(4) .  (3.12) 

k=0 


where  §  is  the  independent  variable  in  Chebyshev  space.  The 
governing  equations  (3.3  and  (3.6),  in  discretized  form,  are  then  arranged 
in  a  generalized  eigenvalue  format.  That  is  if  D  and  E  represent  the 


coefficient  matrices,  then 

DX  =  coEX, 


(3.13) 
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where  the  frequency  o)  =  ac  is  the  eigenvalue,  and  the  eigenvector  X  is 
represented  by 

X  =  [G  H  F  P]T.  (3.14) 

It  should  be  realized  that  E  is  a  singular  matrix.  The  singularity 
is  removed  by  addition  of  the  term  ymp  (which  is  called  an  artificial 
compressibility  factor)  to  the  continuity  equation  (see  Malik  &  Poll 
1985),  where  y  is  a  small  parameter  of  the  order  of  10'^.  The  effect  of 
this  on  the  computed  physical  eigenvalues  is  negligible  as  reported  by 
Khorrami  et  al.  (1989).  The  method  is  global  and  therefore  the  entire 
eigenvalue  spectrum  is  obtained  in  a  single  run.  The  complex  generalized 
eigenvalue  solver  employed  is  the  IMSL  QZ  routine  'EIGZC'. 

The  outer  boundary  conditions  were  enforced  at  rmax  =  100,  and  the 
number  of  Chebyshev  Polynomials  required  varied  depending  on  flow 
conditions,  but  usually  was  in  the  range  between  60  and  80.  At  each 
step,  care  was  taken  to  ensure  that  results  were  at  least  six  or 
seven  significant  figures  accurate.  The  eigenfunctions  were  obtained 
using  an  inverse  Rayleigh's  method  (see  Wilkinson  1965).  The 
discretization  for  the  local  scheme  was  also  spectral;  actually,  the  same 
matrices  D  and  E  were  used  to  compute  the  eigenfunctions. 

Since  it  turns  out  that  results  from  this  study  for  n  =  0  are 
somewhat  different  from  those  of  n  *  0,  this  has  important  implications  on 
the  asymptotic  structure  of  the  solution.  Consequently,  we  shall  consider 
the  axisymmetric  case  separately. 
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4.  Axisvmmetric  (n=0)  modes 

For  n  =  0,  our  numerical  scheme  produced  results  for  wavespeed  eg 
over  a  range  of  values  of  swirl  parameter  q  and  axial  wavenumber  a, 
which  had  the  following  general  features:  (i)  a  number  of  distinct,  real 
modes  exist,  (ii)  all  these  modes  have  eg  <  0,  and  (iii)  these  modes 
were  quite  distinct  from  those  of  other  studies  (e.g.  Lessen  et  al.  1974, 
Duck  &  Foster  1980),  for  which  Cj  *  0.  Indeed  our  routine  was  able 
to  generate  these  other  modes,  which  served  as  a  useful  check  on  the 
accuracy  of  our  scheme.  Results  for  eg  for  the  case  n  =  0,  q  =  1.0, 
over  a  range  of  a  are  presented  in  Fig.l  where  two  distinct  modes  are 
shown.  We  believe  these  to  be  the  two  most  important /dominant  in  this 
case.  We  refer  to  the  mode  represented  by  a  solid  line  as  mode  I,  and  that 
represented  by  a  broken  line  as  mode  II. 

Note  that  the  significance  of  eg  <  0  is  that  no  critical  layers 
exist  (i.e.  <pg  *  0  for  all  r),  a  feature  that  does  lead  to  certain 
simpl i f icat ions . 

The  key  question  now  is  whether  these  modes  are  stable  or  unstable, 
since  the  study  so  far  only  reveals  them  to  be  neutrally  stable  in  the 
limit  of  large  Reynolds  numbers.  To  determine  the  effects  of  viscosity 
on  these  modes,  we  must  consider  terms  OCRg'1)  in  (2.18)  -  (2.21). 

After  some  algebra,  we  obtain  the  following  two  first  order  equations  for 


Gj  and  Pj : 

L1 (Gl ,?1 )  =  Ra  - 

L2(Ci.Pi)  =  Rb- 
We  may  write 

Ra  =  Ra,  +  'Cj  R 


al 


a2  ’ 


(4.1) 

(4.2) 


Rb  =  Rbi  +  iC!  Rbz, 


(4.3) 
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where 


where 


al 


nR3  aRj 

r(P0  <P0 
2W 


Rh  =  -  R2  +  tzt  r3  • 
bi  i  r<pQ  J 


nH0 


a2 


*  •  %  [-r  +  0F°]> 


b2 


r<P0 


Rl  =  -  Fo"  -  }  F0'  +  (  72  +  «2)Fo. 

R2  =  -  GO"  *  7  G0'  +  [  ^  +  “2  ]°0 

.  2n 


H0. 


r3  =  *  h0"  '  7  H0'  +  [  ^  +  “2  ]H0 


^  2n  r 
+  -r  Go- 


(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 


(4.9) 


Here  we  retain  n  since  these  equations  are  also  useful  for  other 
values  of  n. 

In  order  that  (4.1),  (4.2)  have  a  solution,  with  boundary  conditions 
given  by  (2.24),  (2.25),  we  must  have 


1  U°*Ra,  +PtRb,  )dr 


Cl  = 


'O' 


U o+  r.2  +  *  v  )dr 


(4.10) 


>qi  b2 

where  G+  and  P+  are  the  functions  adjoint  to  the  system  (3.8), 
(3.10),  i.e. 


dG+  _  r  n(rW) ’  +  ar2u'  1 

r"  1 — 

-1*0-^]*, 

r3<P0 


(4.11) 
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dP+  _ 
W  ~ 


n2+a2r2 

r2<j>o 


]g+  +  —  P+, 
J  r2<p0 


with  boundary  conditions 


(4.12) 


P+(0)  =  G+’(0)  =  0  if  n  =  0, 

P+'(0)  =  G+(0)  =0  if  n  =  ±1. 

P+(0)  =  G+(0)  =  0  if  ln|>l ,  (4.13) 

and  P+,  G+  -♦  0  as  r  -*  »,  for  all  n. 

Note  that  due  to  the  nature  of  the  solution  for  cq  real,  c\  must  be 
imaginary.  The  adjoint  system  (4 . 1 1 ) - (4 . 12)  is  discretized  similarly  as  in 
the  case  of  the  governing  equations  (3.3)-(3.6).  However,  owing  to  the 
non-linear  occurence  of  <pq  in  (4.11)  which  results  in  a  quadratic  term  for 
the  frequency,  <o,  a  slightly  different  approach  is  taken.  Here,  the 
eigenvalues  were  obtained  using  a  companion  matrix  method.  The  method  is 
very  straight-forward  (see  Bridges  &  Morris  1984  and  Khorrami  et  al. 

1989),  and  involves  linearizing  the  quadratic  term  by  the  following 
transformation: 

P+  -  coP+  =  0.  (4.14) 

which  leads  to  a  third  equation  for  the  adjoint  set.  For  this  case, 
the  eigenvector  X  in  (3.13)  then  becomes 


X  =  [G+  P+  F+]T. 


(4.15) 


It  must  be  mentioned  that  for  each  computation,  the  computed  eigenvalue 
spectrum  of  the  adjoint  system  matched  the  spectrum  associated  with  the 
original  set  i.e.  (3.3)  -  (3.6).  This  is  an  independent  check  on  the 
accuracy  and  integrity  of  our  results. 

To  obtain  C},  a  Gauss  •  Chebyshev  quadrature  is  employed  to 
evaluate  the  integrals  of  (4.10);  the  procedure  is  straightforward. 

To  ensure  accurate  results,  the  number  of  Chebyshev  polynomials,  N, 
was  increased  until  cj  had  converged  to  at  least  five  significant 
figures.  Typically  90  to  100  polynomials  were  more  than  sufficient 
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to  obtain  the  required  accuracy.  Results  for  ©j  =  Im 

for  the  case  n  =  0,  q  =  1.0  are  shown  on  Fig. 2  for  mode  I,  for 

Re  =  5000  and  10,000.  It  is  clear  that  the  first  of  these 

modes  is  destabilised  over  a  range  of  a  with  the  introduction  of  the 

effects  of  viscosity.  Meanwhile,  the  results  for  coj  for  mode  II  are 

shown  in  Fig. 3,  and  it  is  clearly  seen  that  viscosity  stabilizes  this 

mode.  Also  shown  (on  Fig. 4)  are  results  for  coj  =  Im  {acj}  for  mode  I 

obtained  using  the  fully  viscous  routine  of  Khorrami  (1991),  at  the 

same  values  of  q  and  Reynolds  numbers.  We  see  that  the  fully  viscous 

results  appear  to  exhibit  an  upper  neutral  point  which  is  predicted 

extremely  effectively  by  our  asymptotic  results.  The  growth  rate 

acj  in  the  region  of  the  upper  neutral  point  is  also  predicted  accurately 

by  asymptotic  theory.  However,  there  is  an  important  point  of  disagreement 

between  the  Re  »  1  results  and  those  for  the  full  viscous  equations, 

that  concerns  the  nature  of  ocj  as  a  -»  0.  According  to  our  asymptotic 

theory  this  quantity  approaches  a  finite  value  as  a  •*  0,  whilst  the  fully 

viscous  computation  predicts  a  (sharp)  drop  off  at  small  values  of  a. 

However,  this  point  of  disagreement  is  quite  clear.  If  lacjl 
approaches  a  constant  value  as  a  -»  0,  then  C]  =  0(l/a)  as  a  -♦  0, 
and  hence  a  breakdown  in  the  wavespeed  expansion  (3.2)  must  occur. 
Additionally,  if  n  =  0,  it  is  also  clear  that  as  a  -*  0,  <p  -*  0 
(for  bounded  c)  for  all  r.  Specifically  this  breakdown  must 
occur  when  a  =  0(Re_1),  and  hence  we  define  a  scaled  axial  wavenumber 
a  =  Re  a  =  0(1).  (4.16) 

Guided  by  our  previous  results  as  a  -♦  0,  and  by  consideration  of  the 
order  of  magnitude  of  various  terms  in  the  governing  equations, 
for  a  s  0(1)  we  must  have 


F  =  Re  F0(r)  +  0(1), 
G  -  G0(r)  +  0(Re-1). 
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H  =  Re  H0(r)  +  0(1), 

P  =  Re  P0(r)  +  0(1),  (4.17) 

and  the  expansion  for  the  complex  wavespeed  (3.2)  is  retained.  Substituting 
(4.17)  into  (2.18)  -  (2.21)  and  implementing  (4.16),  we  obtain  to  leading 

order 

-  Gq 

Gq'  +  —  +  a  F0  =  0,  (4.18) 

2WH0  . 

-  =  P0r  •  (419) 

-  H0' '  -  I  H0’  +  [  -  idc0  +  idU  +  2^11  ]  H0 

+  [  i  g?  +  i  7  ]  GO  =  0.  (4-20) 

*  F0”  '  7  F0’  +  [  *  *“c0  +  iaU  +  ^  ]  H0 

+  i  G0  +  id  P0  =  0,  (4.21) 

where  the  appropriate  boundary  conditions  may  be  inferred  from 
(2.24)-(2.25).  The  above  system  is  then  an  eigenvalue  problem  for 
CQ(a),  and  is  solved  using  a  simplified  form  of  the  viscous  routine 
used  by  Khorrami  (1991).  However,  owing  to  the  absence  of  the 
eigenvalue  term  in  the  r-momentum  equation,  matrix  E  becomes  singular. 

Here  the  singularity  is  removed  via  a  procedure  which  utilizes  row  and 

column  operations  (see  Metcalfe  &  Orszag  1973).  In  this  procedure, 

the  rank  of  matrices  D  and  E  is  reduced  first  and  the  eigenvalues 

are  then  obtained  using  the  QZ  routine.  Results  for  Im{dcQ(d)} 

are  shown  in  Fig. 5  for  the  case  q  =  1.0,  first  mode.  It  is  quite  clear  that 

a  lower  neutral  point  is  predicted  (at  a  =  90),  whilst  as  a  •*  «, 

Im{aco(d)}  -»  2.3  approximately.  This  is  in  agreement  with  the  values  shown 
on  Fig. 2  as  a  -♦  0.  Indeed  a  routine  asymptotic  analysis  of  the  system 
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(4.18)  -  (4.21)  as  a  -»  ~  confirms  a  correct  asymptotic  match  with 
the  o  =  0(l),  a -»  0  solution. 

Thus,  to  summarise,  we  are  able  to  predict  (using  two  asymptotic 
analyses)  both  the  upper  and  lower  neutral  points  of  this  particular 
unstable  mode  as  well  as  the  temporal  growth  rates  of  the  full 
viscous  equation  results  (Fig. 4).  Furthermore,  the  system 

(4. 1 8)  - (4 . 21 )  was  also  solved  for  the  second  (stable)  mode,  and  it  was 
found  that  this  remained  stable  over  the  entire  range  of  a.  It  should  at 
this  juncture  be  emphasised,  however,  that  as  a  -»  0,  the  use  of  the 
parallel  flow  approximation  is  likely  to  become  increasingly  questionable. 

In  the  following  section  we  go  on  to  consider  the  n  *  0  modes,  paying 
particular  attention  to  cases  for  which  n  =  1,  although  there 
are  some  similarities  with  the  axisymmetric  case,  some  important  and 
interesting  differences  also  exist. 
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S.  Non- ax i symmetric  modes 

The  paper  of  Khorrami  (1991)  presents  results  for  an  unstable  mi  ies 
(with  a  growth  rate  that  diminishes  as  Re  -»  «)  for  the  particular  case 
n  a  1.  In  this  section  we  discuss  the  stability  of  such  an  asymmetric 
disturbance . 

We  initially  follow  the  same  approach  as  that  carried  out  in  the 
previous  sections  and  apply  these  methods  to  the  case  n  =  1, 
q  =  0.7.  Fig. 6  shows  the  variation  of  eg  (which  again  is  real)  with  a, 
obtained  from  the  solution  of  (3.8),  (3.10).  However,  below  a  critical 
value  of  a(  =  ag,  say),  we  see  eg  <  0,  while  above  this  value  eg 
becomes  positive.  We  were  able  to  continue  the  computation  of 
cq  beyond  ag  using  the  numerical  scheme  described  in  Section  3. 
However,  since  cq  >  0,  and  our  numerical  results  suggested  cq  remained 
real,  a  critical  layer  must  be  present  at  any  point  at  which 
90  =  0.  Hence  computation  of  such  modes,  according  to  Lin  (1955), 
must  be  carried  out  by  extending  the  computation  into  complex  r 
space  to  avoid  the  singularity  in  the  differential  equation. 

Although  our  numerical  scheme  performed  extremely  well,  by  its  very  nature 
it  is  not  suitable  for  obtaining  solutions  off  the  real  r  axis.  (The 
authors  did  attempt  a  Runge-Kutta  scheme  for  treating  (3.8),  (3.10), 
but  extremely  small  grid  sizes  which  required  prohibi tat ively  longer 
computer  times,  were  necessary  for  adequate  resolution  of  these  modes, 
even  for  examples  for  which  eg  <  0,  i.e.  for  which  no  critical  layer 
existed.  Further  many  spurious  modes  were  generated  with  this  technique.) 
However,  analysis  below  suggests  why  it  may  be  possible  to  extend  these 
computations  of  (3.8),  (3.10)  into  regimes  where  critical  layers  may 
exist  without  any  special  modification  of  the  scheme,  or  numerical 
di fficul ties. 

Let  us  initially  confine  our  attention  to  values  of  o  <  oq. 
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for  which  the  techniques  and  analysis  of  the  previous  section  are 
applicable  without  modification.  In  particular,  the  computed  values  of 
(Oj  obtained  from  (4.10)  are  shown  on  Fig. 7  as  points  denoted  by  0 
for  the  case  n  =  1,  q  =  0.7,  and  are  to  be  compared  with  the  values 
of  &)j  obtained  using  the  full  viscous  equation  for  the  same  case, 
at  Re  =  10000,  as  presented  in  Fig. 7.  The  computed  and  asymptotic 
values  proved  to  be  indistinguishable  on  the  scale  shown  in  Fig. 7. 

It  is  apparent  that  we  are  able  to  predict,  using  our  asymptotic  theory, 
the  location  of  a  lower  neutral  point  (at  a  =  0.05)  without  the 
requirement  of  a  further  a«  1  substructure.  Agreement  between  the 
asymptotic  and  fully  numerical  results  is  excellent  up  to  a  =  oq, 
the  point  to  which  our  asymptotic  results  extend.  Beyond  a  =  ag, 
the  full  viscous  equations  yield  values  of  gd;  which  are  initially 
seen  to  continue  to  increase  but  then  rapidly  drop  in  value  to  give  an 
upper  neutral  point. 

For  a>  ccq,  it  would  appear  that  there  are  two  distinct  possibilities 
to  explain  this  behaviour. 

The  first,  if  eg  remains  real,  implies  that  a  critical  layer  exists. 
However,  if  we  assume  that  a  is  just  above  the  critical  value  ag.  then 

a  =  ao  +  a.  (5.1) 

where  lal  «  do,  and  o  is  taken  to  be  positive.  As  a  result  we  expect 

co  =  “  cqo  +  0(52) ,  (5.2) 

and  to  be  consistant  with  (3.2),  we  must  have  that  a  »  Re*1.  Suppose 

that  the  location  of  the  critical  layer  is  at  r  =  rg,  then  using  the 

asymptotic  form  for  U(r)  and  W(r),  we  must  have 

r0  *  (  — =*-  )*•  (5.3) 

<*0ac00 

If  this  is  sufficiently  large  then  as  r  -»  rg,  (3.8),  (3.10) 
approximate  to 
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dG0  _  1 
3F  -  r0 


c0 


rg(n2+ag2rg2) 
2n(r-rg)  P°’ 


dPo  2nq  P0 

J_  o0  +  _. 


(5.4) 

(5.5) 


Taking  this  system  (and  indeed  higher  order  terms  that  have  been 

neglected)  reveals  that  the  solution  for  Gq  and  Pg  are  regular 

about  r  =  rg,  i .e. 

00 

Go  =  I  Gon  (r-rg)0.  (5.6) 

n=0 


P()  =  I  P0n  (r-r0)n+1.  (5.7) 

n=0 

and  hence  no  critical  layer  is  required.  This  perhaps  explains  why  we  were 
able  to  extend  our  numerical  scheme  beyond  a  =  oq  without  any  special 
modification  or  difficulties. 

However,  as  a  increases,  rg  moves  toward  the  centre  of  the 
vortex  (r=0).  If  eg  remains  real,  ultimately  the  presence  of  a 
critical  layer  will  become  important^in  particular  its  effect  will 
be  profound  when  rg  =  0(1),  implying  a  -  ag  =  0(1).  Perhaps  it  is  this 
penetration  of  the  critical  layer  close  to  the  vortex  centre  that 
triggers  the  sharp  drop  off  in  growth  rate  acj  with  a,  to  yield 
the  upper  neutral  point  as  seen  in  the  fully  viscous  solutions  in  Fig. 7. 
However,  the  presence  of  the  critical  layer  requires  that  a  detour  be 
made  into  the  complex  r  plane  when  considering  (2.18)  -  (2.21).  According 
to  Lessen  et  al.  1974  this  detour  is  below  the  real  axis  if 
Real  (<p'(rg)}  >  0,  and  vice  versa.  Unfortunately,  as  remarked  earlier, 
our  numerical  scheme  is  confined  to  the  real  axis,  and  so  we  were  unable 
to  carry  out  the  computations  for  our  asymptotic  structure  beyond  a  =  ag 
with  any  degree  of  certainty. 

A  second  possibility  exists,  namely  that  the  inviscid  solution  of  (3.8), 
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(3.10)  yields  complex  stable  values  of  cq  for  a  >  a g.  This 
inviscid  decay  rate  would  then  counteract  the  unstable  viscous 
growth  rate  which  could  result  in  a  rapid  stabilising  trend  for 
Unfortunately,  our  numerical  scheme  was  unable  to  give  a 
categorical  vindication  of  either  of  these  two  possibilities. 


stable 


a  >  a o- 
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6.  Conclusion 

In  this  paper  we  have  presented  asymptotic  analyses  which  describe 
and  indeed  confirm  the  additional  modes  of  instability  recently  found 
numerically  for  the  ' * t rai 1 ing-1 ine  vortex’’  by  Khorrami  (1991). 

These  modes  are  very  different  in  nature  from  those  reported  previously, 
being  inviscidly  neutral,  but  are  clearly  shown  in  this  paper  to 
be  destabilised  by  viscosity.  Previously  reported  modes  of  instability  have 
been  generally  inviscidly  unstable. 

Although  our  investigation  has  been  confined  exclusively  to  the 
' ’ t rai 1 ing- 1 ine  vortex' ' ,  there  is  no  reason  why  such  mechanisms 
should  not  operate  in  the  same  way  for  other  vortex  flows  . 

Finally,  it  must  be  emphasized  again  that  the  parallel  flow 
approximation  has  been  employed  throughout  this  paper,  and  an 
interesting  extension  of  this  work  would  be  to  include  the 
effects  of  non-parallelism. 
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Figure  3.  u>i  vs.  a,  n  =  0,  q  =  1.0,  asymptotic  results,  Re  =  10,000,  mode  II,  —  Numerical, 
- Asymptotic. 


©j  x  104 

2. 


Figure  4. 
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